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Abstract 

Using event driven molecular dynamics simulations, we study a three dimensional one-component 
system of spherical particles interacting via a discontinuous potential combining a repulsive square 
soft core and an attractive square well. In the case of a narrow attractive well, it has been shown 
that this potential has two metastable gas-liquid critical points. Here we systematically investigate 
how the changes of the parameters of this potential affect the phase diagram of the system. We 
find a broad range of potential parameters for which the system has both a gas-liquid critical point 
C\ and a liquid- liquid critical point C2. For the liquid-gas critical point we find that the derivatives 
of the critical temperature and pressure, with respect to the parameters of the potential, have the 
same signs: they are positive for increasing width of the attractive well and negative for increasing 
width and repulsive energy of the soft core. This result resembles the behavior of the liquid- 
gas critical point for standard liquids. In contrast, for the liquid-liquid critical point the critical 
pressure decreases as the critical temperature increases. As a consequence, the liquid-liquid critical 
point exists at positive pressures only in a finite range of parameters. We present a modified van 
der Waals equation which qualitatively reproduces the behavior of both critical points within some 
range of parameters, and give us insight on the mechanisms ruling the dependence of the two 
critical points on the potential's parameters. The soft core potential studied here resembles model 
potentials used for colloids, proteins, and potentials that have been related to liquid metals, raising 
an interesting possibility that a liquid-liquid phase transition may be present in some systems where 
it has not yet been observed. 

PACS numbers: 61.20.Gy, 61.20.Ja, 61.20.Ne, 61.25.Mv, 64.60.Kw, 64.70.Ja, 64.60.My, 65.20,+w 
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I. INTRODUCTION 




The discovery and investigation of liquid-liquid phase transitions in a one- component 
system is of current interest, since recent experiments for phosphorus 0, 0| show a first- 
order phase transition between two stable liquids in the experimentally accessible region of 
the phase diagram. A liquid-liquid phase transition, ending in a critical point was initially 
proposed to exp lain the anomalous behavior of network- forming liquids such as H 2 0,0, 
IIlSlllllIipIilIilIilIMP,IilIi3[il- In particular, the density anomaly, consisting 
in the expansion under isobaric cooling of these systems, has been related to the possible 
existence of a phase transition between low density liquid (LDL) and high density liquid 
(HDL). Simulation results and experimental studies of water predict a LDL- 
transition in an experimentally inaccessible region of the phase diagram || Il2l 
Computer simulations of realistic models of carbon j2j|, phosphorus |22j |. Si02 
(HI, |25[ strongly suggest the existence of first-order LDL-HDL phase transitions in these 
substances. Recently the step changes of the viscosity of liquid metal, such as Co, have been 
theoretically interpreted as evidence of liquid-liquid phase transitions [26[ . 

The presence of the first-order phase transitions in solids and solid-solid critical points, 
determined experimentally j27j and with simulations [H, EE ES EH, have suggested the 
possibility of the existence of liquid-liquid critical points and polymorphism in the amorphous 
state [13, EH Ell • It nas been proposed that systems with solid polymorphism may exhibit 
several liquid phases with local structures similar to the local structures of various crystals. 
Experimental evidence of sharp structural transitions between liquid polymorphs of Se, S, 
Bi, P, I2, Sn, Sb, As2Se3, AS2S3, and M gaBi ? are consistent with phase diagrams with first- 
order liquid-liquid phase transitions j3 a T|3 ^l , analogous to the liquid-liquid phase transition 
seen in rare earth aluminate liquids [36L l37j|. 

These results call for a general interpretation of the basic mechanisms underlying the 
liquid-liquid phase transition. Here we aim to delineate the conditions ruling the accessibility 
of the two liquid phases. A first step in this direction was taken in Refs. [38|, |3£|, where 
we have shown that a specific isotropic soft-core attractive potential, for a one-component 
system, has a phase diagram with LDL-HDL phase transition, with two fluid-fluid critical 
points, and with no density anomaly. 

Here we extend this analysis by varying the parameters of this potential (Fig. [I}. We find 
that, for a wide range of parameters, this potential has a phase diagram with a liquid-liquid 
critical point, and we show how the phase diagram depends on the parameters. We develop 
a modified van der Waals equation (MVDWE) able to describe the behavior of the two 
critical points as a function of the potential parameters, elucidating a mechanism for the 
liquid-liquid phase transition and the conditions under which the liquid-liquid critical point 
occurs at positive pressure. 

In Sec. II we introduce the isotropic soft core potential; in Sec. Ill we describe the 
two different molecular dynamics techniques we use; in Sec. IV we present our results for 
different combinations of parameters that give rise to a liquid-liquid phase transition ending 
in a liquid-liquid critical point; in Sec. V we construct a modified van der Waals equation 
which can qualitatively reproduce the behavior of the two critical points; in Sec. VI we 
discuss the role of potential parameters in changing the position of the critical points; in 
Sec. VII we summarize our results; in Appendix A we present our simulation results for a 
simple square well potential. 



II. THE ISOTROPIC SOFT-CORE ATTRACTIVE POTENTIAL 

For attractive potentials with a sufficiently broad interaction distance, the phase diagram 
has a first-order gas-liquid transition ending in a gas-liquid critical point, and a first-order 
liquid-solid phase transition j4£j. When the attractive range is small, the liquid phase and 
the gas- liquid critical point are metastable with respect to the solid phase (4ll.l42l l43l El. |45| . 
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For a strictly repulsive soft-core potential, simulations show a phase diagram with a 
first-order gas-solid phase transition and a first-order phase transition between two solids of 
different densities, but with the same structural symmetry, ending in a solid-solid critical 
point [H, IH, EH EH- Recent theoretical work has suggested that systems with a broad 
soft-core potential have a fluid-fluid phase transition and liquid anomalies 46], or give rise 
to stripe phases in two dimensions [47]]. 

We have shown in Ref. [38] that the combination of a repulsive soft core with an attractive 
well is sufficient to give rise to a phase diagram with two liquid phases. This simple isotropic 
model potential is similar to those used in the seminal work of Stell and Hemmer j48[, who 
studied a soft-core potential in one dimension (ID). Similar potentials were studied in 2D 
and 3D showing phase diagrams with a possible liquid-liquid critical point [49t l50j] . 

The 3D isotropic potential we consider (Fig. ^) has a hard core (infinite repulsion) at 
distance a, a repulsive soft core of width wr and energy Ur > 0, and an attractive square well 
of width wa and energy —U a < (H,|3^. The potential has three parameters: Wu/a, WA/a, 
and Ur/Ua, where a and U a have been chosen as units of length and energy, respectively. 
Though this potential is discontinuous, it is similar to model potentials for complex fluids, 
such as colloids, protein solutions, star polymers [44], |5l|, |52|, HH H3, HH, |56[ , and resembles 
pair potentials proposed for water [52|, or that have been related to liquid metals under 
specific conditions [53, [H, H^] . 

This potential with parameters wr/ci = 1.0, wa/cl = 0.2, and Ur/Ua = 0.5 has a phase 
diagram with gas-LDL and gas-HDL first-order phase transitions, each ending in a critical 
point in the supercooled fluid region [38[ . Both liquid phases are metastable with respect to 
a single crystal phase and no density anomaly is observed [39j | . 

In this paper we present systematical MD studies of the phase diagrams for this potential 
(Fig. [T]). By varying the parameters of the potential, wa/o,, Wr/cl, and Ur/Ua, we relate 
the attractive and repulsive components of the potential to the appearance and stability of 
the liquid-liquid phase transition and critical points. 



III. MOLECULAR DYNAMICS SIMULATIONS 

We perform MD simulations of N = 850 particles of unit mass m at constant volume 
V, and constant temperature T, interacting via the potential described above (Fig. 1). The 
details of the event-driven MD we use are presented in Refs. [38|, [3£| • We measure time in 

units of all j^^m 1 / 2 and record potential energy and pressure every At = 100 time units. 
To understand the effect of each parameter on the phase diagram of our system, we simulate 
sixteen sets of potential parameters (Table [!}. After a preliminary screening, we choose 
to study the region of parameter space where the low-density, gas-liquid critical point C\ 
always has a critical temperature above that of the high-density critical point Therefore, 
while in Refs. [H, |3£| C2 is a gas-HDL critical point, here C 2 is a LDL-HDL critical point. 
As shown in in Refs. [H, [3£|, C 2 can lie in the supercooled metastable phase, close to the 
line of homogeneous nucleation, as in water or silica 0, EH IH| • We make certain that all 
our calculations are performed before the onset of crystallization, as discussed in Ref. [39]. 
The description of the crystal phases goes beyond the goal of this work. To optimize our 
analysis we use two different MD methods. 



A. Isothermic method 

The first method is a straightforward calculation of the phase diagram's state points. 
For each state point with given p = N/V and T, we perform typically ten independent 
simulations of t ~ 2 x 10 3 time units. We estimate the error in pressure measurements from 
the standard deviation of the ten averaged values computed for each independent simula- 
tion. The state points along the isotherms are approximated by a two-variable polynomial 
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P(p, T) = J2i K a inP %r P K obtained by the least squared fit of all the state points in the vicinity 
of the critical point. This fitting implies mean field critical exponents [6(| and may produce 
incorrect results in the close vicinity of the critical point. However, this method helps us fit 
the state points, known with statistical errors, by approximate polynomial isotherms and 
thus obtain the approximate position of the critical point. 

The coexistence curves are calculated using Maxwell's equal area construction and spin- 
odal line is estimated by locating the maxima and minima of the isotherms. After calculating 
the state points, isotherms, coexistence curves, and spinodal lines, we estimate the critical 
pressure, temperature, and density for G\ and C2 {Pci , Tq\ , pc x , Pc 2 > Tc 2 , and pc 2 , respec- 
tively) as the point where coexistence and spinodal curves meet, coinciding at their maxima. 
We apply this method to six sets of potential parameters (ii, iii, iv, xi, xii and xiv in Tabled). 
The results are presented in Figs. 2-7 in the pressure-density (P-p) phase diagrams. The 
estimates of the critical points are presented in Table |H] 



B. Isochoric method 

The isothermic method gives us fairly complete information about the details of the phase 
diagrams, but requires much computation to calculate enough state points for accurate 
isotherms. Thus, in order to find the positions of critical points for a wide range of potential 
parameters, we adopt a faster but less accurate MD method. For sets of parameters close 
to the sets of parameters studied with the isothermic method, we estimate the location 
of the spinodal line by evaluating the intersections of isochores in the P — T plane. We 
first equilibrate several configurations at a high initial temperature kBTi/UA = 2.0 for 
several values of density above and below the densities where we expect to find pc x and 
pc 2 - At constant density, the system is slowly cooled down from T; to a final temperature 
IcsTf/IIa = 0.1 during a simulation time of 10 4 time units [6]J. 

The average values of T and P are recorded each 100 time units, which is comparable 
to the equilibration time of the system for ksT /Ua > 0.5. As the temperature decreases, 
the equilibration time increases and the method becomes less reliable. Thus, we use this 
method to estimate pressure and potential energy for ksT /Ua > 0.5. 

The error bars of each measurement are of the order of the non-monotonic jumps of the 
isochores (see Fig. El inset). The intersection is determined by fitting isochores with smooth 
polynomial fits. The best results can be achieved by quadratic fits in the temperature range 
including the region of possible isochore crossing extending from 0.9Tq to 1.5Tc, so that the 
tentative critical temperature Tq is inside this interval. 

Since at the spinodal line (dP/dp)x — 0, two isochores with two close values of density 
must intersect in the vicinity of the spinodal line. By definition, the critical point cor- 
responds to the maximum temperature on the spinodal. Therefore, the critical pressure 
and temperature can be evaluated by estimating the pressure corresponding to the max- 
imum temperature at which isochores intersect. The critical density can be estimated as 
(Pi +P2)/2, where p\ and p 2 are the densities of the two isochores intersecting at the highest 
temp eratu re (Fig. [SJ). The critical point values estimated with this method are presented in 
Table Em 

This approximate method is allowed as long as we use it to estimate the critical points 
of potentials with sets of parameters close to those for which we have done a detailed study 
using the isothermic method. We apply the isochoric meth od t o 16 sets of the potential 
parameters. The comparison of the two methods (Tables ITTlllllJ) shows that the resulting 
estimates of critical P, T, and p of C\ and C2 are consistent. 
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IV. PHASE DIAGRAM RESULTS 



Our results in Figs. I2H3 clearly show that the phase diagram strongly depends on the 
potential parameters. For example, phase diagrams in Figs. I2HU have fluid phases (gas, LDL, 
and HDL) at positive pressures, while for the phase diagrams in Figs. EH3 the high-density 
critical point appears at negative pressures, i.e. in the region of stretched fluid. 

To investigate how the position of critical points depends on the potential parameters, 
we vary one of the three parameters wa/cl, wr/cl, Ur/cl at a time, kee ping the other two 
constant. The behavior of T, P, and p for C\ and C2 (Fig. Eland Table llVj) are presented 
in the following. 



A. Effect of the square-well width wa 

By keeping %/a = 0.5 and Ur/Ua = 2.0 constant, we find (Fig. E^-c, Fig. [TUJ) that 
by increasing well width wa, Pci is almost unaffected, while pc 2 decreases, T Cl and T C2 
increase, Pc 1 increases, while Pc 2 decreases. For wa/cl > 0.7 the LDL-HDL critical point C2 
occurs at negative pressures, as in Fig. 03 Hence, C 2 lies in the stretched fluid region and, 
therefore, it is metastable. In order to have a stable LDL-HDL critical point, the attractive 
distance wa/cl must be sufficiently narrow, so that occurs at positive pressures. A too 
narrow well, however, enhances crystallization [3^, El H2, EH, |44], |45[ so that the high- 
density critical point shifts below the line of spontaneous crystallization, becoming difficult 
to observe. Thus the liquid-liquid critical point is observable in our MD simulations only 
for intermediate values of wa/cl. 

B. Effect of the shoulder width wr 

Increasing the width of the repulsive interaction Wr, while keeping wa/ci = 0.7 and 
Ur/Ua = 2.0 constant, we find (Fig. HJl-f, Fig. ITT|) that both pc\ and pc 2 decrease, Tq 1 
decreases, while T C2 increases, and both P Cl and Pq 2 decrease. For wr/cl < 0.4 the dynamics 
of the system in the vicinity of the expected high density critical temperature become too 
slow and the equilibration time becomes too long, with respect to our simulation time, to 
measure the equilibrium state points with sufficient accuracy. Furthermore, as expected for 
decreasing wr, Tq 2 approaches T = (Fig. Efc), suggesting that C2 disappears for wr/cl = 0. 
At wr/cl > 1.0 the system spontaneously crystallizes at high density without showing a 
second critical point Ci- Hence, the width of the shoulder wr/cl must be of an intermediate 
value for C2 to be observed above the lines of spontaneous crystallization and outside the 
region of very slow dynamics, at least for our choice of wa and Ur. 



C. Effect of the shoulder height Ur 

For wr/ci = 0.5 and wa/ci = 0.9, we increase the repulsive energy Ur and find (Fig. ^jg- 
i, Fig. [T2J) that for increasing Ur, pc x decreases, while pc 2 is almost unaffected, both Tq x 
and Tq 2 decrease, Pq x decreases, while Pq 2 rapidly increases. For Ur/Ua < 2.0 the high- 
density phase transition occurs at very low negative pressures and the fluid phases are highly 
metastable. For Ur/Ua > 4.0 the diffusion in the system in the vicinity of the high density 
critical point becomes markedly slow, due to the soft core becoming less penetrable and 
assuming the role of an effective hard core. Therefore, an intermediate repulsive energy is 
needed to observe C 2 in our MD simulations. 
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V. MODIFIED VAN DER WAALS EQUATION 



To rationalize the dependence of the temperature, pressure and density of the two critical 
points on the potential's parameters, we develope a simple mean field theory that gives rise 
to a modified van der Waals equation (MVDWE), 

P = —, 7 - V, (1) 

l- P B(p,T) p > y> 

that has the same form of the standard van der Waals equation (see Appendix A), but with 
an excluded volume B(p,T) depending on the density and temperature of the state point 
and increasing with wr/cl, and with a strength of attraction A that increases with wa/cl and 
decreases with Ur/Ua- It should be pointed out that a different modification of the van der 
Waals equation [£2| also gives rise to the high density critical point. In contrast with our 
work, Ref. 62] is particularly suitable for density dependent potentials since it assumes a 
constant excluded volume B and a density dependent attractive term A(p). 

For a system with a hard core and a soft core, one can assume that the effective ex- 
cluded volume B(p,T) changes with temperature and density [63]. Indeed, at low densi- 
ties and low temperatures, particle cannot penetrate into the soft core so B(p,T) f=a B 2 
where B 2 = 2ii(a + wr) 3 /3 is the excluded volume associated with the soft core. In con- 
trast, for high densities and high temperatures, particles easily penetrate into the soft core 
and B(p,T) w Bi, where B\ = 2ira /3 is the excluded volume associated with the hard 
core. More specifically, B(p,T) must be an analytical function of its parameters such that 
dB(p, T)/dT < 0, dB(p, T)/dp < 0, 



T 

and 



pnB(p,T)=B 1 , (2) 



B 2 , p < 1/B 2 




lim B(p,T) = { XT i/r\ ^ 1 /r ( 3 ) 
r^o VF ' s \ l/p, l/Bi > p > 1/B 2 , v ; 

from which it follows that B(p, T) < 1/p for any p and T > 0. 

Since in any case van der Waals equation can give us only qualitative agreement with 
reality, we can select any model function B(p, T) which satisfies the above conditions. Never 
the less, it is desirable to select B(p,T) in such a way that it will describe the behavior of 
some physical system for which the analytical solution can be found. One dimensional 
system of particles with a pair potential 

r < B { 

U{r) = { U R , B 1 <r<B 2 (4) 
r>B 2 , 

provides such a solution. Applying the Takahashi method [5?1 HH, we obtain the Gibbs 
potential 

G = -keTNi ln(tf jfeflT/BiPi) (5) 

where N\ is the number of particles, T is temperature, Pi is pressure of the one-dimensional 
system and 

= ^ e -PrB x /k B T _ e -PiB 2 /k B Ty-U R /k B T + e -PiB 2 /k B T _ ^ 

Accordingly V\ = dG/dPi and Si = —dG/dT are the volume and entropy of the one- 
dimensional system, and U\ = G — P\V\ + TS\ is the potential energy for the one di- 
mensional system. The fraction of the soft cores f(p,T) penetrated by the particles is 
f(p,T) = Ui(Pi,T) I (NiUr) where Pi must be determined as a function of p from the 
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equation dG/dPi(Pi,T) — V\ = Ni/p. The value = /(p, oo) is the fraction of the 
soft cores penetrated by the particles in the high-temperature limit in which soft cores 
play no role. It can be computed assuming a Poisson distribution of interparticle distances: 
/oo = 1 — e(- Bl "~- B2 "( 1 / p ~" B2 ). The probability that the soft core does not reflect the neighboring 
particle is equal to the fraction of these two quantities ///oo < 1- In this case, the excluded 
volume is equal to B\. In the opposite case with probability 1 — ///oo, the excluded volume 
is equal to B 2 . Hence, the effective excluded volume 



where, 



B(p,T) = f/f oc B 1 + (l-f/f QO )B 2 , 

j ( e -PiBi/k B T _ e -PiB 2 /k B T^ e -U R /k B T 



/oo *(Pi,r)(l - e (Bi-B 2 )/(l/ P -Bi)) 

and Pi must be found from the equation 

1 _ A; B T (B ie - p ^/^ T - B 2 e- p ^l k * T )e- u nl kBT + B^e-***' 1 **' 
p "~P7 + *(Pi,T) 



(7) 
(8) 

(9) 



Figure IT3k illustrates the behavior of B(p,T) for a particular set of parameters. It is 
clear that B(p,T) satisfies all the physical conditions we impose on the effective excluded 
volume. The modified van der Waals equation (JTJ) has two critical points: one for low 
density p << 1/B 2 and another for high density p ~ 1/B 2 , whose positions on the phase 

diagram of the dimensionless variables T = k B T '/Ur, P = B\P/Ur, and p = Bip depend on 
the dimensionless parameters of the MDVWE: B 2 jB\ and A/(UrBi). Figure ITHb shows a 
P — T diagram with two critical points C\, C 2 , for a particular set of parameters, for which 
the position of the critical points are similar to the positions found in our simulations, i.e. 

T C2 <T Cl . 

Now we can relate the parameters of the Eq. ^ to the potential parameters used in our 
simulations. The parameters B\ and B 2 are increasing functions of the hard core diameter 
a and the shoulder width wr, respectively. The parameter Ur has an identical meaning 
in MVDWE and in simulations. The strength of attraction A is an increasing function of 
wa and a decreasing function of Ur. Indeed, according to the formula of the second virial 
coefficient v 2 for our potential, we have [6^| 



v 2 = B 1 + (l- e- UR/kBT ){B 2 - B x ) + (1 - e UA / kBT ) 



27T 

— (a + WR + w A Y - B 2 



(10) 



For large T, it has the form v 2 = B — A/k B T + 0(T~ 2 ) with A = UaVa — UrVr where 
va and vr are positive quantities with the dimension of a volume depending on a, Wr, and 
wa, B = liniT^oo^ and A = lim.j'^ 00 T(B — v 2 ). Hence, in this limit, the virial expansion 
P = k B Tp + k B Tv 2 p 2 + 0(p 3 ) = k B Tp{\ + Bp) — Ap 2 + 0(p 3 ) can be rewritten in the form 
of the van der Waals equation P = k B Tp/ (1 — Bp) — Ap 2 + 0(p 3 ). 

From the equations above we can derive the functional relation for va and v B in the 
limit T — * oo, that are va = (2%/ 3) (a + w B- + wa) 3 — B 2 and Vr = B 2 — B\. By using 
these relations it is possible to see that in general A is an increasing function of wa and a 
decreasing function of Ur. The derivative OA/Owr may have different sign depending on 
other parameters. Although at finite T these relations could be valid only to the leading 
order, it is reasonable to assume that A increases with wa and decreases with Ur at any T. 

However, to simplify our qualitative study of the MVDWE, we assume the parameters 
A, B 2 and Ur are independent. By varying these parameters one at a time and by relating 
B 2 to Wr, and A only to Wa, we found that the MVDWE predicts that the derivatives of 
the low-density critical point values, T Cl , Pc x > Pd , with respect to each of the parameters of 
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MVDWE, have the same sign and this sign is negative for B 2 (wr) and Ur, which increase 
repulsion, and this sign is positive for A(wa), which increases attraction. These results are 
consistent with the MD results for the low-density critical point. For the high-density critical 
point values the MVDWE predicts for some parameters non-monotonic behaviors and we 
find that there are regions of parameters where the critical values as a function of A, B 
and Ur have the same qualitative behaviors as those found in the simulations as a function 
of wa, wr and Ur, respectively (Fig. IT4"|) . Specifically, we observe monotonic behaviors of 
pc 2 (B 2 ), Tc 2 (A), Tc 2 (B 2 ), and Pc 9 {Ur) which qualitatively coincide with the corresponding 
behaviors in simulations fFigs l9~l and 1141 d.b.e. and i). We find non-monotonic behaviors 
of pc 2 (A), pc 2 (Ur), Tc 2 (Ur), Pc 2 (A), and Pc 2 (B 2 ), which qualitatively coincide with the 
corresponding behaviors in simulations in certain range of parameters (Figs 151 and 1141 a. g.h.c. 
and f). 

These observations indicate that the behavior of the critical points in simulation may 
also become non-monotonic in the range of parameters that we do not explore. For exam- 
ple, Tc 2 {Ur) may start to increase for large Ur/Ua > 4 and small wr/cl < 0.5. Another 
interesting prediction of the MVDWE is that for large B 2 /Bi > B t (ABi/Ur), where Bf{x) 
increases from Bt{0.7) = 1 to 5^(3.2) = 1.7 the high-density critical temperature becomes 
larger than the low-density critical temperature as in simulations of Refs. [38L l39| , for which 
the repulsive shoulder wr/cl = 1 was much wider than the attractive well wa/cl = 0.2. Also, 
MVDWE predicts the existence of the third, very high density critical point for large BjjB^ 
and large ABi/Ur, which was recently observed in simulations with a wide soft core |66l |. 




VI. ROLE OF POTENTIAL PARAMETERS 

In the following we will present the comparison between the MD results and MVDWE 
predictions. 

A. The low density critical point 

First we note that at low densities, corresponding to the critical point C\, and at suffi- 
ciently low temperatures, particles do not penetrate into the repulsive region, r < a + wr. 
Therefore, we can assume that, at low enough temperatures and densities, the system is 
interacting via an effective potential given by a simple square- well with hard core a + wr, 
an attractive well of relative width uu/(a + Wr) and attractive energy U a- 

Indeed, for increasing width of the attractive well Wa, Pc x is roughly constant (Fig. EK) 
and Td and Pc x increase (Fig. and c). This behavior is consistent with the predictions of 
the standard van der Waals theory f or the g as-liquid critical point for a square-well potential 
(see Appendix A), that yields Eqs. (jA7HA9|) . This result supports the idea that the effect of 
the soft core is negligible at low densities. The MVDWE also predicts strong increase of Pc x 
and Td with the strength of attraction A, which increases with wa- For pc 1) the MVDWE 
predicts a weak increase, which can be observed in Fig. for larger wa- 

For increasing width of the repulsive shoulder Wr, pc x decreases and saturates (Fig. 01) 
and Td and Pc x decrease (Fig. and f). As a consequence of the above considerations, the 
increase of wr corresponds to a decrease of this effective attractive parameter. Accordingly, 
Tc\ and Pc\ display the behavior predicted for the decreas e of the attractive width in van 
der Waals theory for the square- well potential (Eqs. IA8HA9|) . Moreover, the behavior of pc x 
is consistent with Eq. (|A7|) . which predicts pc x ~ 1/a 3 . Indeed, in our case, the hard core 
is replaced by the effective hard core, hence a 3 pc\ ~ a 3 / (a + wr) 3 , which is a decreasing 
function of wr. The calculations using MVDWE completely confirms these predictions by 
showing that Tc\,Pc\, and pc x all decrease with increasing wr (or B 2 = 2n(a + wr) 3 /3 as 
in Sec.V). 
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For increasing repulsive energy Ur, the behaviors of pc x , and Pc x are the same as 
those observed for increasing wr (Fig.[S^-i). This can be understood by considering that the 
increase of Ur effectively decreases the penetrability of interparticle distances r < a + wr. 
The soft core becomes an effective hard core as described above. In particular, the saturation 
of pd , already observed in Fig. Eli, is now more evident and an analogous behavior is now 
also seen for Tc 1 and Pc\- This result shows that for a high-enough repulsive energy Ur 
and low-enough T, the soft-core potential is equivalent to a hard-core potential, for which 
there is no dependence of the critical point on Ur. Again, the MVDWE agrees with these 
predictions by showing that T Cl ,Pd, and pc x decrease with increasing Ur. 

B. The high density critical point 

For increasing wa, the critical point density pc 2 decreases, the temperature Tc 2 increases, 
and the pressure Pq 2 decreases (Fig. EJi-c). This finding is in agreement with MVDWE 
predictions for a wide range of parameters. fFig. I14b-c). The behavior of Tc 2 is consistent 
with the idea that the increase of the attractive distance increases the overall attractive 
strength of the potential, allowing more particles to fit within the attractive interaction 
range. As a consequence, the system enters the low-energy and high-density fluid phase at 
a higher temperature, i.e. Tq 2 increases. Hence, the increase of wa increases the average 
kinetic energy of particles at C%, favoring the overcoming of the soft-core shoulder at low 
pressures, and we can expect that P C2 decreases. Moreover, the increase of wa decreases 
the number of elastic interparticle collisions at the soft-core distance, hence decreases their 
contribution to the virial expression for the pressure ( see Eq. (11) in Ref.[39]), decreasing 
the critical pressure Pa, ■ Note that the behavior of Pq 2 in this case is the opposite of the 
behavior for P Cl (Fig. Et). With the decrease of pressure, the density must also decrease, so 
we can conclude that pc 2 must decrease with increasing wa, in agreement with our simulation 
results. 

For increasing wr, T C2 increases and saturates, and both P C2 and pc 2 decrease with a 
tendency toward saturation (Fig. EJi-f), in agreement with predictions of MVDWE for a 
wide range of parameters (Fig. IT4"H-f). This happens because the transition from the LDL 
to the HDL is characterized by the penetration of particles into the repulsive soft cores of 
their neighbors. Therefore the repulsive soft-core distance a + Wr characterizes the typical 
distance between the particles at Ci- Hence the increase of wr reduces the critical density 
Pc 2 - The behavior of pressure follows the behavior of density, as in the case of wa, while 
the derivatives of pressure and temperature must have the opposite sign due to the same 
arguments as above. 

For increasing Ur, Pq 2 increases, pc 2 slowly increases, and Tq 2 decreases (Fig. Efe-i). 
The predictions of MVDWE coincide with the behavior of Pq 2 and pc 2 in a wide range of 
parameters (Fig. IT4"g.i) . However, the theory apparently predicts an increasing Tq 2 with 
Ur, except for very small UrB\/A < 0.4 and large B2/B1 > 1.5 (Fig.ll4h). This discrepancy 
arises from the fact that, although it is physically clear that the attractive strength A is a 
decreasing function of Ur, we find the explicit dependence of A on Ur only in the limiting 
case T — > 00, while for finite T we assume them to be independent. Hence, we ignore that 
an increase of Ur decreases A which induces, as shown in Sec.V, a decrease in Tc 2 . 

The behavior of Pq 2 is easier to understand. Indeed, the pressure at which the repulsive 
shoulder can be overcome increases with Ur, which is consistent with the increase of Pq 2 
with Ur. This effect is expected to be more evident at high Ur and to saturate for decreasing 
Ur, which is consistent with our results. The critical density pc 2 increases with Ur, for small 
values of Ur/Ua, as a consequence of the increase of Pc 2 , and is practically independent of 
Ur when the soft core plays the role of an effective hard core, i. e. for large enough Ur/Ua- 
The decrease of Tc 2 with the increase of Ur is more difficult to explain. Nevertheless, the 
same argument as in the case of wa and wr which predicts that the derivatives Tc 2 and Pc 2 
must have the opposite signs may apply in this case as well. Finally, we note that increasing 
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with Ur and decreasing with wr, the behavior of Pc 2 in three dimensions is consistent with 
its behavior in the one dimensional case, for which P C2 /Ua = (Ur/Ua + ^)/wr [68]. 

VII. DISCUSSION AND CONCLUSIONS 

We have studied an isotropic attractive soft-core square potential in three dimensions 
that has a phase diagram with a gas-liquid critical point C\ and a liquid-liquid critical point 
C2, separating a HDL and a LDL phases. We have investigated, with molecular dynamics 
simulations, how the critical density, temperature and pressure of the two critical points 
vary as a function of the three parameters of the potential, that are the repulsive energy 
Ur/Ua in units of the attractive energy Ua, the repulsive width wr/ci in units of the hard 
core a, and the attractive width w^/a in units of a. 

Table IIVI and Fig. |§] show our results for the p, T, and P of C\ and C2 for varying 
parameters of the potential. To summarize, the behavior of C\ is consistent with that of a 
system interacting via an effective square-well potential, with a hard core a + wr, a relative 
attractive well wa/(cl + wr) and attractive energy Ua- The increase of Ur/Ua or of wr/cl 
decreases the effective attractive strength and this effect saturates for large values of Ur/Ua- 

This behavior is perfectly predicted by the simple mean field MVDWE. In MVDWE, as 
in the standard van der Waals equation, the relevant physical parameters are the excluded 
volume B and the stength of attraction A, but now we assume that B = B(p,T) increases 
with wr/ci and depends on the state point, while A increases with wa/cl and decreases with 
Ur/Ua- 

The MVDWE predictions are consistent also with our results on C2. In general, this 
approach rationalizes why the increase of Ur/Ua and the decrease of wa/ci have the same 
qualitative effect on the critical points, since they have the same effect of the attractive 
strength. Moreover, it rationalizes the effect of increasing wr via the increase of the excluded 
volume, hence the decrease of the critical densities. This decrease induces a decrease of the 
critical pressures Pq x and Pc 2 , as a consequence of the mechanical stability of the fluid 
phases. 

While for the low-density critical point C\, the decrease of Pc 1 occurs with the decrease 
of Tc 1 - I i.e. their derivative with respect to the potential parameters always have the same 
sign, for the high-density critical point C2 the critical pressure and temperatures always 
have derivatives with the opposite sign. This behavior can be understood in terms of the 
number of elastic collisions with the soft core, which decreases as Tc 2 increases, reducing 
the virial contribution to the pressure. At the same time, the increase of Tq 2 reduces the 
pressure Pc 2 necessary to overcome the repulsive soft core and enter into the HDL phase. 

As a consequence, the high density critical point C2 exists at positive pressure only in a 
finite region in the parameter space. Indeed, when the attraction is too strong, i.e. wa/(i 
is too large or Ur/Ua is too small, the pressure Pq 2 becomes negative. On the other hand, 
when the strength of attraction is too week, C2 occurs in the deeply supercooled liquid phase, 
becoming difficult to observe as in the experimental situation of water or silica El HIl • 

In conclusion, the behavior of both low-density and high-density critical points quali- 
tatively obeys the mean field predictions of the modified van der Waals theory based on 
effective excluded volume which varies between the hard-core value for high temperature 
and low density and the soft-core value for low temperature and low density. The quanti- 
tative theory based on the thermodynamic perturbation approximations or various integral 
equation closures |k| is yet to be developed. 

Confirming the results presented in Refs. |38l . l39| | we do not find density anomaly. Our 
simulations show that density anomaly is unlikely to exist for the discontinuous double-step 
potentials shown in Fig. ^ in contrast to ramp potentials [49] and a Gaussian soft core 
potential [531 ]. 

Our results may be relevant for experiments on systems that can be described by an 
isotropic soft-core attractive potential and have no density anomaly, such as colloids, protein 
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solutions or liquid metals. Indeed, our results show that in these systems the possibility of 
the existence a liquid-liquid phase transition will depend on the relative ratio between the 
attractive and the repulsive part. 



Acknowledgments 

We thank C. A. Angell, V. V. Brazhkin, A. Geiger, I. Kohl, T. Keyes, T. Loerting, V. 
Ryzhov, R. Pastor, and F. W. Starr for helpful discussions and the NSF Chemistry Division 
(CHE-0096892) for support. 



APPENDIX A: SQUARE- WELL FLUID SYSTEM 

Here we recall the case of a square well attractive potential, where the only parameter 
is the width of the attractive well wa/o>, since the hard core distance a and the attractive 
energy Ua can be taken as units of distance and energy, respectively. In particular, we 
show how the gas-liquid critical point density p c , temperature T c , and pressure P c depend 
on wa/ci- 

Even in this simple case, the phase diagram has no exact a naly tical solution and one 
must rely on various approximations and numerical simulations [4y, [70, El, EH • Using 
MD simulations of N = 850 particles, we verify that the be havi or of p c , T c , and P c , are 
approximately linear for a wide range of wa/cl (Fig. EE Table [40]. The values of o?p c 
decreases for increasing wa/o, and ci 3 P c /Ua and ksTc/UA increase with wa/cl. 

Except for density, these results are in agreement with the van der Waals theory (see, 
e.g., Ref. |65[). The equation of state in the van der Waals theory is given by 

p^-V, (AD 

where B = |a 3 7r has the meaning of excluded volume per particle and A is a quantity, with 
the dimension of the product of energy and volume, characterizing the strength of attraction 
between particles. Therefore, A can be related to the product of Ua and the volume of the 
attractive well, which is proportional to wao 2 for small wa/o>- 
The position of a critical point must satisfy the equations 



dP 
dp 

d 2 p 



dp 2 



, (A2) 



. (A3) 



For the van der Waals equation of state Eq. (|A1|) . Eqs. (|A2|) and flA3|) yield the coordinates 
of the critical point: 

ft = Jb • (A4 » 

<^ = ||, (A5) 

1 A 
27 'W 
Hence: 

p c o? ~ const , (A7) 



p c=^ir 2 - (A6) 
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k B T c w a 

a3Pc Wa (A9) 



U A a 

which predict that k B T c /UA and o?P c /Ua increase with wa/cl, while p c a 3 does not depend 
on it. 
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TABLE I: Sets of parameters for the generic soft-core potential (Fig. [IJ considered in this paper: 
wr/ci and wa/cl are the soft-core width and the attractive width, respectively, both in units of the 
hard-core distance, and Ur/Ua is the repulsive energy in units of the attractive energy. Sets i-vi 
have same wa and Ur; sets ii, vii-xii have same wr and Ur; sets xii-xvi have same wr and wa- 



Set a 


UUr/(1 


wa/cl 


Ur/Ua 


i 


0.4 


0.7 


2 


ii* 


0.5 


0.7 


2 


hi* 


0.6 


0.7 


2 


iv* 


0.7 


0.7 


2 


V 


0.8 


0.7 


2 


vi 


0.9 


0.7 


2 


vii 


0.5 


0.3 


2 


Vlll 


5 


4 


2 


ix 


0.5 


0.5 


2 


X 


0.5 


0.6 


2 


xi* 


0.5 


0.8 


2 


xii* 


0.5 


0.9 


2 


xiii 


0.5 


0.9 


2.5 


xiv* 


0.5 


0.9 


3 


XV 


0.5 


0.9 


3.5 


xvi 


0.5 


0.9 


4 


a The asterisk (*) 


denotes sets for which critical points 


are calculated via two methods (see 


Tables HH and 



ED. 



TABLE II: Temperatures Tq 1 and Tq 2 , pressures Pc x and Pc 2 ; and densities pc x and pc 2 > for the 
critical points C\ and C2, respectively, computed by the isothermic method. 



Set 


k B T Cl /U A 


a 3 P Cl /U A 


a 3 Pd 


k B T C2 /U A 


a 3 P C2 /U A 


a PC 2 


ii 


1.30 ±0.01 


0.04 ±0.01 


0.11 ±0.02 


0.58 ±0.02 


0.15 ±0.02 


0.33 ±0.02 


iii 


1.24 ±0.01 


0.03 ±0.01 


0.09 ± 0.02 


0.69 ± 0.02 


0.11 ±0.02 


0.28 ±0.02 


iv 


1.18 ±0.03 


0.025 ± 0.003 


0.08 ± 0.02 


0.75 ±0.01 


0.07 ±0.01 


0.24 ±0.02 


xi 


1.52 ±0.01 


0.05 ±0.01 


0.11 ±0.02 


0.69 ±0.01 


-0.11 ±0.01 


0.33 ± 0.02 


xii 


1.82 ±0.01 


0.06 ± 0.02 


0.12 ±0.02 


0.96 ±0.02 


-0.21 ±0.02 


0.32 ±0.03 


xiv 


1.59 ±0.01 


0.043 ± 0.004 


0.10 ±0.02 


0.58 ±0.01 


-0.01 ±0.01 


0.35 ±0.02 
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TABLE III: Temperatures Tq x and Tq 2 , pressures Pq 1 and Pq 2 , and densities pc x and pc 2 for the 
critical points C\ and C2, respectively, estimated by cooling the system at constant p (isochoric 
method) for the potential with the set of parameters in Table [I] 



Set 


k B T Cl /U A 


a 3 P Cl /U A 




a 3 Pd 


k B T C2 /U A 


a 3 P C2 /U A 


a 3 pc 2 


i 


1.34 ±0.02 


0.04 ±0.01 





13 ±0.02 


0.47 ±0.01 


0.28 ±0.01 


0.42 ± 0.03 


ii 


1.32 ±0.01 


0.04 ±0.02 





11 ±0.02 


0.62 ± 0.02 


0.19 ±0.02 


0.33 ± 0.02 


iii 


1.25 ±0.01 


0.03 ±0.01 





09 ±0.01 


0.69 ± 0.02 


0.11 ±0.01 


0.29 ±0.02 


iv 


1.19 ±0.01 


0.03 ±0.01 





08 ±0.01 


0.74 ± 0.01 


0.07 ±0.01 


0.26 ±0.02 


V 


1.15 ±0.02 


0.02 ±0.02 





07 ±0.01 


0.75 ±0.01 


0.04 ±0.01 


0.22 ±0.02 


vi 


1.11 ±0.02 


0.02 ±0.02 





07 ±0.01 


0.76 ±0.01 


0.03 ±0.01 


0.20 ±0.02 


vii 


0.68 ± 0.01 


0.02 ±0.01 





12 ±0.01 


0.48 ± 0.03 


2.22 ± 0.02 


0.46 ± 0.06 


viii 


0.82 ±0.01 


0.03 ±0.01 





12 ±0.01 


0.52 ± 0.03 


1.65 ±0.02 


0.42 ± 0.03 


ix 


0.96 ±0.01 


0.03 ±0.02 





11 ±0.01 


0.53 ± 0.03 


1.05 ±0.03 


0.39 ± 0.05 


X 


1.12 ±0.01 


0.04 ±0.01 





10 ±0.01 


0.57 ±0.01 


0.58 ±0.01 


0.35 ±0.02 


xi 


1.54 ±0.02 


0.05 ±0.02 





12 ±0.01 


0.70 ±0.01 


-0.09 ±0.01 


0.33 ±0.03 


xii 


1.84 ±0.02 


0.06 ±0.02 





13 ±0.01 


0.96 ±0.01 


-0.22 ±0.01 


0.31 ±0.03 


xiii 


1.67 ±0.01 


0.05 ±0.01 





11 ±0.01 


0.72 ±0.01 


-0.15 ±0.01 


0.35 ±0.01 


xiv 


1.62 ±0.02 


0.05 ±0.01 





09 ±0.01 


0.60 ± 0.01 


0.01 ±0.01 


0.37 ±0.04 


XV 


1.57 ±0.01 


0.04 ±0.01 





09 ±0.01 


0.55 ±0.01 


0.28 ±0.01 


0.35 ±0.02 


xvi 


1.54 ±0.01 


0.04 ±0.01 





09 ± 0.01 


0.53 ± 0.02 


0.60 ± 0.02 


0.35 ±0.02 



TABLE IV: Summary of the effects on pc±, Pc x , and pc 2 , Tc 2 , Pc 2 , from variation of 

parameters w A /a, wr/cl, and Ur, one at the time. The symbols f, J. and « represent, respectively, 
an increase, a decrease and a small variation of a thermodynamic quantity as a consequence of the 
increase of the potential parameter. 
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TABLE V: Density p Cl temperature T c , and pressure P c for the gas-liquid critical point for a single 
square- well potential with attractive range w A , attractive energy U A and hard-core distance a. 



w A /a a 3 p c k B T c /U A a 3 P c /U A 

0.1 0.45 ±0.05 0.48 ±0.01 0.062 ± 0.005 

0.2 0.42 ±0.05 0.67 ±0.01 0.076 ± 0.005 

0.3 0.34 ±0.05 0.86 ±0.01 0.088 ± 0.005 

0.4 0.35 ±0.05 1.04 ±0.01 0.094 ± 005 

0.5 0.30 ±0.03 1.24 ±0.01 0.103 ±003 

0.6 0.28 ±0.03 1.45 ±0.01 0.113 ±003 



TABLE VI: Single square-well potential: parameters of the linear fits in Fig. ^] for k B T c /U A = 
p + qw A / 'a, and analogous linear expressions for a 3 P c /U A and a 3 p c , for the gas-liquid critical point. 





k B T/U A 


a 3 P/U A 


a 3 p 


V 


0.29 ±0.01 


0.055 ± 0.002 


0.482 ± 0.008 


(1 


1.91 ±0.03 


0.097 ± 0.006 


-0.35 ±0.02 
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FIG. 1: The generic soft-core potential with attractive well with parameters wa/ci, wr/o-, and 
Ur/Ua- We use the parameters listed in Table [I] 
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FIG. 2: The MD P-p phase diagram for the potential in the inset, with the parameter set ii in 
Tabled The long-dashed lines are the fits of the calculated state points (circles) at constant T. The 
isotherms (from top to bottom) are for k B T/U A = 1-30, 1.29, 1.28, 1.27, 1.26, 1.25, 1.24 at low-p and 
ksT/UA = 0.62, 0.60, 0.58, 0.57, 0.55, 0.53, 0.50 at high-p. The fits are calculated by considering 
P a polynomial function of both T and p. The isotherms show two regions with negative slope, 
i.e. mechanically unstable, delimited by the spinodal lines (solid bold lines). Each spinodal line 
is associated with a first-order phase transition. By using the Maxwell construction, we estimate 
the coexisting regions associated to each spinodal line, delimited by the phase transition line (bold 
dashed line). The coexisting regions are clearly separated at the considered temperatures. The 
phase transition line at low p is indistinguishable from the spinodal line at this scale. The points 
where the coexisting lines merge with the spinodal lines are, by definition, the critical points C\ (at 
low-p) and C2 (at high-p). No spontaneous crystal nucleation is observed in the explored region of 
the phase diagram. 
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FIG. 3: As in Fig. [21 for parameter set iii in Table [IJ The isotherms in the low-p region (from 
top to bottom) are for k B T/U A = 1-25, 1.24, 1.23, 1.22, 1.20, 1.18, 1.16, and in the high-p region 
are for kBT/UA = 0.72, 0.70, 0.68, 0.65, 0.62, 0.60. Spontaneous crystal nucleation is observed for 
T < T C2 and p> pc 2 - 
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0.1 0.2 0.3 

Density a 3 p 

FIG. 4: As in Fig. |2l for parameter set iv in Table d The isotherms (from top to bottom) in 
the low-/? region are for k^T jXJ a = 1-20, 1.15, 1.10, 1.05, 1.00, and in the high-/? region are for 
ksT/UA = 0.77, 0.75, 0.73, 0.72, 0.70, 0.69. No spontaneous crystal nucleation is observed in the 
explored region of the phase diagram. 
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FIG. 5: As in Fig. [21 for parameters xi in Table [I] The isotherms (from top to bottom) in the 
low-/? region are for ksT/UA = 1.53, 1.52, 1.515, 1.51, 1.50, 1.48, 1.46, and in the high-p region 
are for k B T/U A = 0.70, 0.68, 0.66, 0.64, 0.63, 0.62, 0.61. C 2 is at negative pressure. 
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0.1 0.2 0.3 0.4 

Density a 3 p 

FIG. 6: As in Fig. [2 for parameters xii in Table HJ The isotherms (from top to bottom) in the 
low-/? region are for k B T/U A = 1-83, 1.82, 1.815, 1.81, 1.80, 1.79, 1.75, 1.70, and in the high-p 
region are for kBT/UA = 0.98, 0.96, 0.64, 0.92, 0.90. C2 is at a negative pressure. 
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FIG. 7: As in Fig. [21 for parameters xiv in Table HJ The isotherms (from top to bottom) in the 
low-/? region are for ksT/UA = 1.65, 1.62, 1.60, 1.58, 1.55, 1.50, 1.45, and in the high-p region are 
for k B T/U A = 0.60, 0.59, 0.58, 0.57, 0.56, 0.54. C 2 is at negative pressures. 
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FIG. 8: Estimation of the critical point C2 by the isochoric method for the set of potential pa- 
rameters wa/cl = 0.5, wr/ci = 0.5, Ur/Ua = 2. Inset: P at constant a? p = 0.492 for the MD 
calculation during the slow cooling described in the text. For kBT/UA > 0.5 the errors on the 
estimate of the state points are of the order of the non-monotonic jumps. The interpolating line 
is a quadratic fit of the calculated points, and gives an estimate of the isochore at a 3 p = 0.492 
for k B T/U A > 0.5. Main panel: quadratic fits of isochores for a 3 p = 0.492, 0.435, 0.405, 0.387, 
0.361 (from top to bottom). The critical point C2 is located at the highest-T intersection of two 
isochores (region inside the circle). The indeterminacy of this intersection gives an estimate of the 
error on the values T C2 = 0.53 ± 0.03, Pc 2 = 1.05 ± 0.03, and pc 2 = 0.39 ± 0.05. 
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FIG. 9: The behavior of the density, temperature and pressure of the low-density critical point 
C\ (open circles) and high-density critical point C2 (filled squares) for variations of the potential 
parameters (a)-(c) wa, (d)-(f) w R and (g)-(i) U R . The other two parameters are constant: in 
(a)-(c) w R /a = 0.5 and U R /U A = 2, in (d)-(f) w A /a = 0.7 and U R /U A = 2, in (g)-(i) w A /a = 0.9 
and w R /a = 0.5. Where not shown, errors are smaller than the symbol size. Lines are guides for 
the eye. 
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FIG. 10: The gas-LDL critical point {C\) and LDL-HDL critical point (C2) in the P — T plane, for 
varying attractive width wa and constant wr/ci = 0.5, Ur/Ua = 2.0. Symbols denote: wa/cl = 0.3 
(circles), 0.4 (left triangles), 0.5 (diamonds), 0.6 (up triangles), 0.7 (right triangles), 0.8 (down 
triangles), 0.9 (squares). Open symbols are for C\ and filled symbols are for C2. The arrows 
denote the direction of increasing wa- 
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FIG. 11: The gas-LDL critical point (Ci) and LDL-HDL critical point (C2) in the P — T plane, 
for varying shoulder width wr/ci and constant wa/cl = 0.7, Ur/Ua = 2.0. Symbols denote: 
wr/ci = 0.4 (circles), 0.5 (up triangles), 0.6 (diamonds), 0.7 (left triangles), 0.8 (down triangles), 
0.9 (squares). Open symbols are for C\ and filled symbols are for C2. The arrows denote the 
direction of increasing wr. 
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FIG. 12: The gas-LDL critical point (Ci) and LDL-HDL critical point (C2) in the P — T plane, 
for varying repulsive energy Ur/Ua and constant wa/cl = 0.9, wr/cl = 0.5. Symbols denote: 
Ur/Ua = 2.0 (circles), 2.5 (up triangles), 3.0 (diamonds), 3.5 (down triangles), 4.0 (squares). Open 
symbols are for C\ and filled symbols are for C2. The arrows denote the direction of increasing 
Ur. 




FIG. 13: (a) The behavior of the effective excluded volume B(p,T) for a one dimensional system 
with .Ba/Bi = 2 for densities B lP = 0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9 from top to bottom. The 
thick curve indicates the behavior for B\p = 0.5 = B1/B2. (b) The isochores of MVDWE for 
B2/B1 = 1.4, A/{B\Ur) = 2.2 on the P — T plane. An open circle indicates the low density 
critical point. A filled circle indicates high-density critical point. A dashed line indicates the 
low-density critical isochore B\pc 1 ~ 0.25. A dot-dashed line indicates the high-density critical 
isochore B\pc 2 ~ 0.70. The thick line indicates B\p ~ 0.71 = B\jBi- Note that isochores start to 
develop density anomaly below the high-density critical point. 
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FIG. 14: The behavior of the density, temperature and pressure of the high-density critical point 
C2 for variations of the MVDWE parameters (a)-(c) A, (d)-(f) B 2 and (g)-(i) Ur. The other two 
parameters are constant: in (a)-(c) B2/B1 = 1.1 (o), 1.5 (□), 2.0 (o), 3.0 (A) and Ur/Ua = 1-0, in 
(d)-(f) A/iBtfjR) = 0.2 (o), 1.0 (□), 2.0 (o), 3.0 (A) and U R /U A = 1, in (g)-(i) A/(BiUr) = 1.0 
and B 2 /Bi = 1.1 (o), 1.3 (*), 1.5 (□), 2.0 (o), 3.0 (A). Lines are guides for the eye. 



28 




FIG. 15: Inset: the single square-well potential denned by the well width wa- Main panel: symbols 
represent the values of the normalized temperature T c (circles), the density p c (diamonds), and 
the pressure P c (squares) of the gas- liquid critical point for different values of the parameter wa- 
Where not shown, errors are smaller than the symbol size. Lines are the linear fits T c , P c , and p c 
as functions of wa, with the parameters in Table [VTl 
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